#################

# Causes of perceived government trustworthiness
# European Journal of Political Research
# Replication code -- please see the readme.

#################

# set up --------------------------------------------------------------------

remove(list=ls())

install.packages("pacman")
library(pacman)

p_load(here, rio, ggpubr, dplyr, jtools, estimatr, grf, MASS, stargazer,
       reporttools, nnet, margins, panelr, cregg, lmtest, FindIt, viridis)

all_conjointdf <- import(here("data", "all_conjointdf.RData"))

# ---- Figure 3 -----

mm_all <- cj(all_conjointdf, chosen ~ performance + interests + transparency +
               tenure + ideology + priority + approval, estimate = "mm", 
             h0 = 0.5, id = ~new_id,
             feature_labels = list("performance" = "Competence", 
                                   "interests" = "Benevolence",
                                   "transparency" = "Integrity",
                                   "tenure" = "Tenure",
                                   "ideology" = "Ideology", 
                                   "priority" = "Priorities",
                                   "approval" = "Approval"))

mm_all %>% 
  ggplot(aes(level)) +
  geom_hline(yintercept = 0.5, color="grey10") +
  geom_linerange(aes(y=estimate, ymin= lower, ymax = upper, color = feature),
                 linewidth=3, position = position_dodge(width = 0.9), alpha = 0.7) +
  geom_point(aes(y = estimate, color = feature), position = position_dodge(width=0.1), size = 4) + 
  facet_grid(feature~., scales='free', space='free') +
  coord_flip() +
  labs(x=NULL, y='Marginal means') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.spacing = unit(2,'pt'),
    panel.background = element_rect(color='black'),
    legend.position = "none",
    legend.title = element_blank(),
    strip.background= element_blank(),
    strip.text = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.text.x = element_text(size = 16),
    axis.text.y = element_text(face="bold", size = 10)) +
  guides(fill = guide_legend(reverse=T)) +
  scale_fill_viridis(discrete = TRUE, option = "plasma", 
                     begin = 0.15,
                     end = 0.8, direction = 1) +
  scale_colour_viridis(discrete = TRUE, option = "plasma", 
                       begin = 0.15,
                       end = 0.8, direction = 1) +  
  geom_text(
    aes(label = sprintf("%0.2f", round(estimate, digits = 2)), x = level, y = 0.35, color = feature),
    position = position_dodge(width = 0.9), size = 4, fontface = "bold",
  )

ggsave(here("figures", "all_mm.png"),
       width = 21, height = 24, units = "cm",
       dpi = 600)

# ---- Figure 4 -----

m1_mm <- cj(all_conjointdf, chosen ~ performance + interests + transparency +
              tenure + ideology + priority + approval, estimate = "mm", 
            h0 = 0.5, by = ~country, id = ~new_id,
            feature_labels = list("performance" = "Competence", 
                                  "interests" = "Benevolence",
                                  "transparency" = "Integrity",
                                  "tenure" = "Tenure",
                                  "ideology" = "Ideology", 
                                  "priority" = "Priorities",
                                  "approval" = "Approval"))
m1_mm %>% 
  ggplot(aes(level)) +
  geom_hline(yintercept = 0.5, color="grey10") +
  geom_linerange(aes(y=estimate, ymin= lower, ymax = upper, color = feature),
                 linewidth=3, position = position_dodge(width = 0.9), alpha = 0.7) +
  geom_point(aes(y = estimate, color = feature), position = position_dodge(width=0.1), size = 4) + 
  facet_grid(feature~country, scales='free') +
  coord_flip() +
  labs(x=NULL, y='Marginal means') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.spacing = unit(2,'pt'),
    panel.background = element_rect(color='black'),
    legend.position = "none",
    legend.title = element_blank(),
    strip.background= element_blank(),
    strip.text = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.text.x = element_text(size = 16),
    axis.text.y = element_text(face="bold", size = 10)) +
  guides(fill = guide_legend(reverse=T)) +
  scale_fill_viridis(discrete = TRUE, option = "plasma", 
                     begin = 0.15,
                     end = 0.8, direction = 1) +
  scale_colour_viridis(discrete = TRUE, option = "plasma", 
                       begin = 0.15,
                       end = 0.8, direction = 1) +
  ylim(0.3, 0.65)

ggsave(here("figures", "country_mm.png"),
       width = 32, height = 25, units = "cm",
       dpi = 600)

# ---- Figure 5 -----

ame <- CausalANOVA(formula = chosen ~ performance + interests + transparency,
                    data = all_conjointdf, 
                    pair.id = all_conjointdf$new_id,
                    diff = T,
                    cluster = all_conjointdf$task,
                    ord.fac = c("TRUE", "TRUE", "TRUE"),
                    nway = 2,
                    screen=T)

ame <- summary(ame)[[4]]
limit <- max(abs(ame$AMIE)) * c(-1, 1)

ggplot(ame, aes(y= Level2, x=Level1)) +
  geom_raster(aes(fill=AMIE)) +
  theme_classic() +
  scale_fill_gradientn(colours=c("#5601A4FF", "white", "#FCA636FF"),
                       breaks=c(-0.04, -0.05, -0.03, -0.02, -0.01, 0, 0.01, 0.02, 0.03, 0.04, 0.05), limit=limit) +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  labs(x = "First level",
       y = "Second level")

ggsave(here("figures", "amies.png"),
       dpi = 600)

# ---- Figure 6 ---- 
# NB: can take a while to run 

devtools::install_github("tsrobinson/cjbart")

# keep only the core variables, outcome, and id. 

bart_df <- all_conjointdf %>% dplyr::select(-c(choice, profile, task, id))


cj_model <- cjbart::cjbart(data = bart_df,
                           Y = "chosen", 
                           id = "new_id") 

het_effects <- cjbart::IMCE(data = bart_df, 
                            model = cj_model, 
                            attribs = c("performance","interests","transparency","tenure","ideology", "priority", "approval"),
                            ref_levels = c("Competent","Usually acts in interests of citizens",
                                           "Open and transparent","Elected very recently","Center", 
                                           "Economy", "30% approve, 70% disapprove"),
                            cores = 1)

hetvimp_m1 <- cjbart::het_vimp(het_effects)

res <- hetvimp_m1$results

res <- res %>% mutate(attribute = case_when(Level %in% c("Competent", "More competent than incompetent",
                                                           "More incompetent than competent", "Incompetent") ~ "Competence",
                                              Level %in% c("Usually acts in interests of citizens", "Interests of majority, not of some",
                                                           "Some citizens, but not all", "Its own interests") ~ "Benevolence",
                                              Level %in% c("Open and transparent", "Usually open and transparent",
                                                           "Usually not open and transparent","Not open and transparent") ~ "Integrity",
                                              Level %in% c("Elected very recently","Elected fairly recently","Elected a long time ago") ~ "Tenure",
                                              Level %in% c("Education", "Environment",  "Immigration",  "Welfare", "Economy","Health") ~ "Priority",
                                              Level %in% c("Very left wing", "Fairly left wing","Center-left","Center","Center-right",
                                                           "Fairly right wing", "Very right wing") ~ "Ideology",
                                              Level %in%  c("30% approve, 70% disapprove","45% approve, 55% disapprove","60% approve, 40% disapprove") ~ "Approval"))

limit <- max(abs(res$importance)) * c(-1, 1)

ggplot(res, aes(y= Level, x=covar)) +
  geom_raster(aes(fill=importance)) +
  facet_wrap(~attribute, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme_minimal() +
  scale_fill_gradient(low = "white", high = "red",
                      breaks=c(0, 75, 125)) +
  theme(axis.text.x = element_text(hjust=1, angle = 45),
        strip.text.x = element_text(size = 9)) +
  labs(x = "Respondent level",
       y = "Attribute level") +
  guides(fill=guide_legend(title="Importance")) +
  scale_x_discrete(labels=c("Age", "Country", "Education", "Employment", "Sex", "Ideology", "Trust"))

ggsave(here("figures", "interactions.png"), 
       width = 21, height = 22, units = "cm",
       dpi = 600)


# ---- Appendix ----

task_check <- cj(all_conjointdf, chosen ~ performance + interests + transparency +
                   tenure + ideology + priority + approval, estimate = "mm", 
                 h0 = 0.5, by = ~ task, id = ~new_id,
                 feature_labels = list("performance" = "Competence", 
                                       "interests" = "Benevolence",
                                       "transparency" = "Integrity",
                                       "tenure" = "Tenure",
                                       "ideology" = "Ideology", 
                                       "priority" = "Priorities",
                                       "approval" = "Approval"))


plot(task_check, group = "task", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5)

ggsave(here("figures", "appendix", "task_check1.png"), 
       width = 21, height = 22, units = "cm")

taskcheck2 <- all_conjointdf %>% dplyr::filter(country == "Argentina") %>% droplevels(all_conjointdf$task) %>%
  cj(., chosen ~ performance + interests + transparency +
       tenure + ideology + priority + approval, estimate = "mm", 
     h0 = 0.5, by = ~task, id = ~new_id,
     feature_labels = list("performance" = "Competence", 
                           "interests" = "Benevolence",
                           "transparency" = "Integrity",
                           "tenure" = "Tenure",
                           "ideology" = "Ideology", 
                           "priority" = "Priorities",
                           "approval" = "Approval"))

task_arg <- plot(taskcheck2, group = "task", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="Task", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5)

taskcheck3 <- all_conjointdf %>% dplyr::filter(country == "Croatia") %>% droplevels(all_conjointdf$task) %>%
  cj(., chosen ~ performance + interests + transparency +
       tenure + ideology + priority + approval, estimate = "mm", 
     h0 = 0.5, by = ~task, id = ~new_id,
     feature_labels = list("performance" = "Competence", 
                           "interests" = "Benevolence",
                           "transparency" = "Integrity",
                           "tenure" = "Tenure",
                           "ideology" = "Ideology", 
                           "priority" = "Priorities",
                           "approval" = "Approval"))

task_cro <- plot(taskcheck3, group = "task", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="Task", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5)

taskcheck4 <- all_conjointdf %>% dplyr::filter(country == "Spain") %>% droplevels(all_conjointdf$task) %>%
  cj(., chosen ~ performance + interests + transparency +
       tenure + ideology + priority + approval, estimate = "mm", 
     h0 = 0.5, by = ~task, id = ~new_id,
     feature_labels = list("performance" = "Competence", 
                           "interests" = "Benevolence",
                           "transparency" = "Integrity",
                           "tenure" = "Tenure",
                           "ideology" = "Ideology", 
                           "priority" = "Priorities",
                           "approval" = "Approval"))

task_esp <- plot(taskcheck4, group = "task", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="Task", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5)

taskcheck5 <- all_conjointdf %>% dplyr::filter(country == "France") %>%
  cj(., chosen ~ performance + interests + transparency +
       tenure + ideology + priority + approval, estimate = "mm", 
     h0 = 0.5, by = ~task, id = ~new_id,
     feature_labels = list("performance" = "Competence", 
                           "interests" = "Benevolence",
                           "transparency" = "Integrity",
                           "tenure" = "Tenure",
                           "ideology" = "Ideology", 
                           "priority" = "Priorities",
                           "approval" = "Approval"))

task_fr <- plot(taskcheck5, group = "task", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="Task", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5)

ggarrange(task_fr + 
            theme(axis.text.y = element_blank(),
                  axis.ticks.y = element_blank(),
                  axis.title.y = element_blank(),
                  strip.text = element_blank() ), 
          task_arg + theme(axis.text.y = element_blank(),
                           axis.ticks.y = element_blank(),
                           axis.title.y = element_blank(),
                           strip.text = element_blank(),
                           legend.position = "none"), 
          task_cro + theme(axis.text.y = element_blank(),
                           axis.ticks.y = element_blank(),
                           axis.title.y = element_blank(),
                           strip.text = element_blank(),
                           legend.position = "none"),
          task_esp+ theme(axis.text.y = element_blank(),
                          axis.ticks.y = element_blank(),
                          axis.title.y = element_blank(),
                          strip.text = element_blank(),
                          legend.position = "bottom"),
          common.legend = T,
          ncol = 2,
          nrow = 2)

ggsave(here("figures", "appendix", "task_check_all.png"), 
       width = 21, height = 22, units = "cm")

all_conjointdf <- all_conjointdf %>% mutate(
  profile = as.factor(profile)
)

prof_check <- cj(all_conjointdf, chosen ~ performance + interests + transparency +
                   tenure + ideology + priority + approval, estimate = "mm", 
                 h0 = 0.5, by = ~profile, id = ~new_id,
                 feature_labels = list("performance" = "Competence", 
                                       "interests" = "Benevolence",
                                       "transparency" = "Integrity",
                                       "tenure" = "Tenure",
                                       "ideology" = "Ideology", 
                                       "priority" = "Priorities",
                                       "approval" = "Approval"))


plot(prof_check, group = "profile", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5)

ggsave(here("figures", "appendix", "prof_check1.png"), 
       width = 21, height = 22, units = "cm")

profcheck2 <- all_conjointdf %>% dplyr::filter(country == "Argentina") %>% droplevels(all_conjointdf$task) %>%
  cj(., chosen ~ performance + interests + transparency +
       tenure + ideology + priority + approval, estimate = "mm", 
     h0 = 0.5, by = ~profile, id = ~new_id,
     feature_labels = list("performance" = "Competence", 
                           "interests" = "Benevolence",
                           "transparency" = "Integrity",
                           "tenure" = "Tenure",
                           "ideology" = "Ideology", 
                           "priority" = "Priorities",
                           "approval" = "Approval"))

prof_arg <- plot(profcheck2, group = "profile", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="Profile", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5) +
  labs(title = "Argentina")

profcheck3 <- all_conjointdf %>% dplyr::filter(country == "Croatia") %>% droplevels(all_conjointdf$task) %>%
  cj(., chosen ~ performance + interests + transparency +
       tenure + ideology + priority + approval, estimate = "mm", 
     h0 = 0.5, by = ~profile, id = ~new_id,
     feature_labels = list("performance" = "Competence", 
                           "interests" = "Benevolence",
                           "transparency" = "Integrity",
                           "tenure" = "Tenure",
                           "ideology" = "Ideology", 
                           "priority" = "Priorities",
                           "approval" = "Approval"))

prof_cro <- plot(profcheck3, group = "profile", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="Profile", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5) +
  labs(title = "Croatia")

profcheck4 <- all_conjointdf %>% dplyr::filter(country == "Spain") %>% droplevels(all_conjointdf$task) %>%
  cj(., chosen ~ performance + interests + transparency +
       tenure + ideology + priority + approval, estimate = "mm", 
     h0 = 0.5, by = ~profile, id = ~new_id,
     feature_labels = list("performance" = "Competence", 
                           "interests" = "Benevolence",
                           "transparency" = "Integrity",
                           "tenure" = "Tenure",
                           "ideology" = "Ideology", 
                           "priority" = "Priorities",
                           "approval" = "Approval"))

prof_esp <- plot(profcheck4, group = "profile", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="Profile", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5) +
  labs(title = "Spain")

profcheck5 <- all_conjointdf %>% dplyr::filter(country == "France") %>%
  cj(., chosen ~ performance + interests + transparency +
       tenure + ideology + priority + approval, estimate = "mm_differences", 
     h0 = 0.5, by = ~profile, id = ~new_id,
     feature_labels = list("performance" = "Competence", 
                           "interests" = "Benevolence",
                           "transparency" = "Integrity",
                           "tenure" = "Tenure",
                           "ideology" = "Ideology", 
                           "priority" = "Priorities",
                           "approval" = "Approval"))

prof_fr <- plot(profcheck5, group = "profile", size = 4.6, feature_headers = F)+
  facet_wrap(~feature, ncol = 1L, scales = "free_y", strip.position = "right") +
  theme(legend.position = "right") +
  guides(colour=guide_legend(title="Profile", reverse = T)) +
  theme(strip.background= element_blank(),
        strip.text = element_text(size = 17),
        axis.text.y = element_text(face="bold", size = 17),
        legend.text=element_text(size=12)) +
  geom_vline(xintercept = 0.5, linetype="dotted", alpha = 0.5) +
  labs(title = "France")

ggarrange(prof_fr + 
            theme(axis.text.y = element_blank(),
                  axis.ticks.y = element_blank(),
                  axis.title.y = element_blank(),
                  strip.text = element_blank() ), 
          prof_arg + theme(axis.text.y = element_blank(),
                           axis.ticks.y = element_blank(),
                           axis.title.y = element_blank(),
                           strip.text = element_blank(),
                           legend.position = "none"), 
          prof_cro + theme(axis.text.y = element_blank(),
                           axis.ticks.y = element_blank(),
                           axis.title.y = element_blank(),
                           strip.text = element_blank(),
                           legend.position = "none"),
          prof_esp+ theme(axis.text.y = element_blank(),
                          axis.ticks.y = element_blank(),
                          axis.title.y = element_blank(),
                          strip.text = element_blank(),
                          legend.position = "bottom"),
          common.legend = T,
          hjust = -0.8,
          #labels = c("France", "Argentina", "Croatia", "Spain"),
          ncol = 2,
          nrow = 2)

ggsave(here("figures", "appendix", "prof_check_all.png"), 
       width = 21, height = 22, units = "cm")

#### Tables #####

# Pooled MMs

esttab <- dplyr::select(mm_all,
                        "Feature" = feature,
                        "Level" = level,
                        "Est." = estimate,
                        "SE" = std.error,
                        "LowerCI" = lower,
                        "UpperCI" = upper)

export_mm <- bind_cols(esttab) %>%
  arrange(as.numeric(Feature), -as.numeric(Level)) %>%
  dplyr::select(Level, Est., SE, LowerCI, UpperCI)

xtable::xtable(export_mm, 
               caption = "Marginal means", 
               digits = 3) -> esttab

xtable::print.xtable(esttab, 
                     include.rownames = F,
                     tabular.environment = "longtable",
                     type = "latex",
                     path = here("tables"),
                     file = "mm_table2.tex")


# Country-separated MMs
esttab <- dplyr::select(m1_mm,
                        "Country" = country,
                        "Feature" = feature,
                        "Level" = level,
                        "Est." = estimate,
                        "SE" = std.error,
                        "LowerCI" = lower,
                        "UpperCI" = upper)

export_mm <- bind_cols(esttab) %>%
  arrange(as.numeric(Feature), -as.numeric(Level)) %>%
  dplyr::select(Country, Level, Est., SE, LowerCI, UpperCI)

xtable::xtable(export_mm, 
               caption = "Marginal means", 
               digits = 3) -> esttab

xtable::print.xtable(esttab, 
                     include.rownames = F,
                     tabular.environment = "longtable",
                     type = "latex",
                     #path = here("tables"),
                     file = "mm_table.tex")

# BES MMs

esttab <- dplyr::select(uk_mm,
                        "Feature" = feature,
                        "Level" = level,
                        "Est." = estimate,
                        "SE" = std.error,
                        "LowerCI" = lower,
                        "UpperCI" = upper)

export_mm <- bind_cols(esttab) %>%
  arrange(as.numeric(Feature), -as.numeric(Level)) %>%
  dplyr::select(Level, Est., SE, LowerCI, UpperCI)

xtable::xtable(export_mm, 
               caption = "Marginal means from Figure \ref{fig:ukmm}", 
               digits = 3) -> esttab

xtable::print.xtable(esttab, 
                     include.rownames = F,
                     tabular.environment = "longtable",
                     type = "latex",
                     path = here("tables"),
                     file = "ukmm_table.tex")


